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CONVERSION FACTORS, U.S. CUSTOMARY TO METRIC (SI) UNITS OF MEASUREMENT 


U.S. customary units of measurement used in this report can be converted to 
metric (SI) units as follows: 


Multiply by To obtain 
inches 25.4 millimeters 
Dy BYP centimeters 
square inches 6.452 square centimeters 
cubic inches 16.39 cubic centimeters 
feet 30.48 centimeters 
0.3048 meters 
square feet 0.0929 Square meters 
cubic feet 0.0283 cubic meters 
yards 0.9144 meters 
Square yards 0.836 Square meters 
cubic yards 0.7646 cubic meters 
miles 1.6093 kilometers 
square miles 259.0 hectares 
knots 1.852 kilometers per hour 
acres 0.4047 hectares 
foot-pounds 1.3558 newton meters 
ainlabares T0197 x 10723 kilograms per square centimeter 
ounces 28.35 grams 
pounds 453.6 grams 
0.4536 kilograms 
ton, long 1.0160 metric tons 
ton, short 0.9072 metric tons 
degrees (angle) 0.01745 radians 
Fahrenheit degrees 5/9 Celsius degrees or Kelvins! 


1To obtain Celsius (C) temperature readings from Fahrenheit (F) readings, 
use formula: C = (5/9) (F -32). 
To obtain Kelvin (K) readings, use formula: K = (5/9) (F -32) + 273.15. 


SYMBOLS AND DEFINITIONS 


Fourier series coefficients 

distance from bottom to pressure sensors 

wave celerity 

cospectrum value 

total water depth 

breaking wave depth 

wave energy density 

complex Fourier coefficient 

discrete frequency value 

ratio of rms breaking wave height to breaking wave depth 
acceleration of gravity 

breaking wave height 

counting indexes 

dynamic pressure response factor 

wave number 

wavelength 

sensor spacing 

index to account for gage number 

total number of discrete data points 

frequency number, argument of Fourier series coefficients 
longshore energy flux at the surf line 

pressure time-series values 

dynamic pressure 

quad—-spectrum value 

ratio of unwindowed energy density to windowed energy density 
complex cross-spectrum value 

length of time series record 

high frequency cutoff period 

weighting coefficient 

water surface elevation 

gage orientation angle 

specific weight of seawater 

wave direction 

average mean depth of water overlaying pressure sensors 
frequency step 

time step 


angular wave frequency 


COMPUTER ALGORITHM TO CALCULATE LONGSHOKE ENERGY FLUX AND 
WAVE DIRECTION FROM A TWO PRESSURE SENSOR ARRAY 


by 
Todd L. Walton, Jr. and Robert G. Dean 


I. INTRODUCTION 


The documented (FORTRAN IV programing language) computer program discussed 
in this report was originally written as part of the Coastal Engineering 
Research Center's (CERC) Longshore Sand Transport Research Program and was used 
in analysis of wave data collected at Channel Islands Harbor in conjunction 
with a study of sand transport at Channel Islands Harbor as discussed in Bruno, 
et al. (1981). 


The program performs the basic analysis of two wave gage pressure records 
necessary to compute wave direction and wave energy at a given frequency and 
computes the longshore energy flux used in sand transport for the entire energy 
spectrum of the wave record. This program uses linear wave theory for the wave 
transformation process and includes the assumption of straight and parallel 
bottom contours necessary for application of Snell's law of refraction. 


Necessary steps in the analysis of the wave data are presented in Sections 
II and III of this report. Subroutines are discussed and sample outputs for 
some wave records from the Channel Islands wave gage pressure sensor pair are 
given. 


The program presently accepts data in the standard CERC magnetic-tape for- 
mat where record lengths consist of 4,100 values. The first four values are 
the gage number and the date-time group, and the remaining 4,096 values are 
the pressures recorded in thousandths of a foot (head) of water at 0.25-second 
intervals. Should other input data be available, the program could easily 
be modified to accept the data by simple changes in the main program and in 
subroutines BUF and SWITCH. 


Sample outputs have been presented for real wave data; some wave direc— 
tional information cannot be obtained for all frequencies because the spectral 
information at some frequencies is ill-conditioned. The percent of energy for 
which this problem occurs is a small part of the energy (usually <3 percent) of 
the entire spectrum and is insignificant in energy-flux computations. Reasons 
for this feature are discussed later. 


Il. METHODOLOGY 
Calculating the longshore energy flux at breaking required the following 
steps: 


(1) Calculation of the frequency—by-frequency wave direction and 
energy at the location of the wave gages; 


(2) determination of the breaking wave depth; 


(3) transformation of the wave spectrum to the “breaker” line, 
including shoaling and refraction effects; and 


(4) computation of ecm the longshore energy flux at the 
surfline. 


Each of the steps is described below. 


1. Calculation of Wave Direction and Energy Spectrum at Wave Gages. 


As noted previously, each of the input time-series pressure records con- 
sists of 4,096 data points with a time increment of 0.25 second. To reduce 
computational costs, modified time series are formed for analysis by averaging 
four adjacent data points. These new time series contain 1,024 data points 
spaced at 1.0-second intervals. This increases the aliasing period from 0.5 
to 2.0 seconds; however, this is justified as the pressure response factor 
for a water depth of 6 meters and a wave period of 2 seconds is approximately 
0.005. 


The time series are analyzed using a standard fast Fourier transform (FFT) 
program to determine the coefficients. For example, for pressure time series 
from gage 1 


oe i2tmnj 
P\Gj) = J fay(n) ~ ab, (n)] exp (2771) eo) 


n=0 


in which i =/-1 and N is the total number of data points, T/At = 1,024, 
where T is the time series record length of 1,024 seconds, At the time 
increment of 1 second between samples, and j a discrete time t, where tj = 
discrete time value = jAt. The FFT coefficients are defined in terms of the 
pressure time series as 


N- 


1 
aj(n) = iby(n) =F J Py(3) exp (-1 272d) (2) 

j=0 
where the argument "n"” of the Fourier coefficients a(n) and b(n) speci- 
fies the quantity to be a discrete function of wave frequency, fa where 
f,» a discrete frequency value, is nAf (where Af = 1/T) and the a,(0) term 
represents the mean value of the time-series pressure record for wave gage l. 
Similar relationships exist for wave gage 2. In calculating the FFT coeffi- 
cients, there are several options that may be employed in an attempt to reduce 
spectral leakage which arises due to representing an aperiodic time series by 
a periodic series. A large number of possible data windows (weighting func- 
tions for data) have been developed to reduce the adverse effects of spectral 
leakage (Harris, 1974). These can be expressed in the form of a weighting 


function w(j), such that the modified time series p'(j) is of the form 
p'(j) = w(j) pi) 


in which p(j) is the digitized measured pressure value at time tj = jAt, and 
w(j) a weighting function. A characteristic of these weighting functions is 
that they are equal to unity at the midpoint of the time series and decrease 
to a lesser value near the two ends. In the present program, a “cosine bell” 
weighting function is used; however, through comparisons of P g With and 
without this function, it was established that the effect of the weighting 
function was minimal (<5 percent). The cosine bell weighting function is 
expressed by 


w(j) = 1 (1.0 - cos 7) (3) 


It is clear that the application of a weighting function will reduce the total 
energy in the record. This effect is partly compensated for by the following 


equation: 
"G) pe "() (4) 
p'(j) = /——— PNG 
<p'2> 


thereby ensuring the same total energy in the altered and original time 
series, where <p> is the mean square value of the original time series 
and <p'2> the mean square value of the weighted time series. It is the 
altered time series p"(j) that is subjected to FFT analysis. The primes 
will be dropped hereafter for convenience. The average mean depth of water 
overlying the pressure sensors, Ad, is obtained by averaging the m time 
series to obtain a, (0) - For two separate time series records, m = 1, 2 (wave 
gages 1 and 2), 


Ad = 0.5 [a (0) + a9(0)] (5) 


The total water depth, d, is the sum of Ad and the distance, B, of 
the pressure sensors above the bottom (in later examples B = 0.76 meter). 


Each FFT pressure coefficient is transformed to a water surface displace- 
ment coefficient by the following linear wave theory relationship discussed in 
the Shore Protection Manual (SPM) (see Ch. 2, U.S. Army, Corps of Engineers, 
Coastal Engineering Research Center, 1977): 


water surface coefficients dynamic pressure coefficients 


1 
[a (), bn (2) 1, = Ga [a,(), b,(2) 1, (6) 


Zz 
in which the subscripts n and p denote water surface and dynamic pressure 
coefficients, respectively. The factor 


ae cosh k(n) B (7) 
a ~ aan Tm) al 


where y is the specific weight of fluid (seawater) and is included when 
pressure coefficients are in normal units of pressure (i.e., N/M2 or equiva-— 
lent). Im equation (7), B represents the distance of the pressure sensors 
above the bottom and k(n) is the wave number associated with the angular 
frequency, w(n) = (2mnAf), as obtained from the linear wave theory dispersion 
relationship 


w(n)* = gk(n) tanh k(n) d (8) 


One of the disadvantages of measuring waves with near—bottom pressure sen- 
sors is evident by examining equations (6) and (7). For the higher frequen- 
cies (shorter wave periods) K,(n) is very small which means that the higher 
frequency waves result in very small pressure fluctuations near the sea floor. 
Thus, to avoid contaminating the calculated water surface displacements, it is 


usually necessary to apply a high frequency cutoff, above which the pressure 
contributions are discarded. The proper selection of this high frequency 
cutoff depends on the signal to noise characteristics of the pressure sensor 
and the signal conditioning system. In the present program, the high fre- 
quency cutoff was established at a wave period of 3.0 seconds. Wave gage 
analyses by Thompson (1980) have shown that a 3.0-second high frequency 
spectral cutoff value provides reasonable estimates of total wave energy at 
west coast (U.S.) locations. 


Denoting hereafter the FFT coefficients for the water surface as a(n) 
and b(n), it is noted that the coefficients have the following properties: 


N- 
<n@> = , [a2(n) + b*(n)] (9) 
and ee 
a($+n)-a(Z-n) (10) 
b(Z+ n)= -b(4- 2) (11) 
and thus N/2 
Gie> SPD ) fiaeG) 2 be G@)] (12) 


n=] 


Thus, the total (kinetic and potential) energy E(n) associated with a par- 
ticular wave frequency component, n, is 


E(n) = 2y[a*(n) + b2(n)] (13) 


Now consider two wave or pressure sensors located at (x) > y,) and (x > Yond) 
(see Fig. 1). The results will be developed considering discrete frequencies. 


= Bet 
B Beta Shoreline 


( parallel to shoreline ) 


Gage 2 (x2, ye) yo 


(normal to bottom contours )} 


Gage t (x,,y;) 


Figure 1. Definition sketch for two sensor array. 
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The water surface displacement consistent with the assumption of one 
direction per frequency is 


N-1 
n(x, ys j) = ) Fm) exp filnw,t - k(n) x - ky(m) y]} 
: ee 
(14) 
Nol i2tnj 
= ) [a(@) = ib@)] exo ( ) 
n=0 N 


where 0, is the primary analysis frequency (= 2m/record length = 21/T 
2mA£), and O(n) the direction of wave propagation at frequency w(n) = nw, + 
The wave number components, k, (2) and (n), are expressed in terms of the 


wave number, k(n), and wave direction, O(n), as 


k, (n) = k(n) cos O(n) (15) 


k(n) = k(n) sin Q(n) (16) 


The cross spectrum, $1542)» of the two measured water surface displacements 
(or dynamic pressures) is given by 


Si2(m) = |F(n) |2 {exp - i [k(n) cos O(n)(x2 - x)) 
(17) 
+ k(n) sin O(n) (y, = yds 


Denoting the separation distance and angle as 2 and 8, respectively, the 
cross spectrum can be expressed as (see Fig. 1) 


S}2(n) |F(n) |? {cos [k(n) cos (O(n) - 8)] - i sin [k(n) &cos (O(n) - B)]} 


cospectrum (n) — i quad-spectrum (n) (18) 


Cj2(n) - iQi2(n) 


Thus, from equation (18), the wave direction O(n) associated with each wave 
frequency can be expressed as 


1 Q,,() 


@(n) = B F cos ! / ————. tan — 
k(n) & C1 ,(n) 


(19) 


The above relationship has two roots, one of which must be selected based 
on physical considerations of the most likely direction of wave propagation. 
In the present case, assuming no wave reflection from the beach, the ambiguity 
in wave direction is ruled out; for wave sensors nearly parallel to the beach, 
the minus sign in equation (19) is appropriate. 
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There are two conditions for which it was not possible to calculate the 
wave directions O(n). These include poorly conditioned wave data, presumably 
due to spectral leakage, and spatial aliasing due to large separation distance 
between the two gages. If the data are poorly conditioned for determining 
wave direction, the absolute value of the quantity within the brackets {-} in 
equation (19) may exceed unity, a physically impossible condition since the 
extreme values of the cosine function are +l. This tends to occur for the 
extremely long waves for which the energy is small and the value of k(n) is 
also small, the latter tending to result in large values of the bracketed 
quantity. The percentage of energy for which this condition occurred in the 
analysis of one year's wave data collected at Channel Islands Harbor was 
relatively small, averaging 2 to 3 percent with a maximum of approximately 
10 percent. The second condition is related to spatial aliasing and requires 
that one-half the wavelength be equal to or greater than the projection of the 
Wave gage separation distance in the direction of wave propagation. Referring 
to Figure l, 


L > 22 {cos[@(n) - B]} (20) 


max 
which indicates that for the least adverse effects of spatial aliasing, the 
gages should be on an alinement parallel to the dominant orientation of the 
wave crests. As will be discussed later, in calculating Peg an attempt was 
made to account for this effect of aliasing by augmenting the calculated 
values, illustrated as follows by 


Eror 
C2) ean = Ug are (21) 


in which the subscripts c and cm indicate calculated and calculated modi- 
fied, respectively. E and E represent the total wave energy values 
and the wave energy not affected by spatial aliasing or poorly conditioned 
wave data, respectively. The total wave energy is that energy in the wave 
spectrum below the high frequency spectral cutoff value. 


2. Transformation of Wave Spectrum to Breaker Line. 


At this stage, the wave energy and wave direction in the vicinity of the 
gages are determined. These values are then transformed to the breaker line 
accounting for wave refraction and shoaling. 


To determine the wave breaking depth, the onshore-directed energy flux is 
calculated in accordance with the expression (based on Snell's law of refrac- 
tion) and equated to an equivalent expressed in terms of wave characteristics 
at breaking. 


N/2 
Onshore energy flux = ) y2 [a(n)* + b(n)2] C,(n) cos O(n) 
i= 
(22) 
_ YEE 
= mer Cop cos O 
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Assuming that the breaking wave angle, O,, is small, that the waves will 
break under shallow-water conditions, and that the ratio of breaking wave 
height to depth is a constant, the breaking wave height, Hp» is then given 
by 


—(N/2 \ CB \o*2 
Hy ={ ), 16 la@)e = b(n) 2] %(n)iFeos!o(a))\02* = (23) 
n=] : J 8 


where GB is the ratio of root-mean-square (rms) breaking wave height to 
breaking depth, GB = Hy /dy (here assumed GB = 0.78). With the breaking depth 
known, each wave component is transformed to shore accounting for both wave 
refraction and shoaling based on linear wave theory. 


Wave refraction is in accordance with Snell's law and the assumption 
that straight and parallel contours existed between the gage and breaking 
locations 


C_.(n) 


-1 C, 4 
0,(n) = sin sin 0,.(n) (24) 
ie 


where C is linear wave celerity (see the SPM, Ch. 2) in which the r_ sub- 
scripts denote the “reference (gage) location. 


With the wave energy and direction now known at the breaker line, the 
value of the longshore energy flux, (Pp.).,, is readily determined 


(Pes dom = R(Peg)e 


(25) 
N/2 


R}2y ) [a*(n) + b7(n)], [C,(n)], [cos O(n) sin O(n)], 
n=] 


in which the factor R is given by the ratio 


Eror 
ny Ss —— 
E 


as defined in and discussed in relation to equation (21). 


III. MAIN PROGRAM DOCUMENTATION 


The detailed programing steps in analysis for the longshore energy flux, 
CPs em? (which in this program is calculated in terms of rms wave height) are 
presented in this section. Program steps are numbered to correspond to areas 
in the program listing where computations are carried out. A program listing 
with corresponding numbered steps follows the program documentation. Note 
that preceding text has used the indexes j and n for time and frequency, 
respectively, while the program which follows uses the index I for both time 
and frequency. A listing of the main program is presented in Figure 2. 
Program steps are as follows and refer to numbered parts of main program 
listing: 


13 


10 


20 


2s 


30 


38 


40 


4§ 


50 


$3 


60 


65 


70 


anaonr 
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PRUGKAM SPECTC INPUT» OUTPUT» TAPES@INPUT» TAPEGSOUTPUTs TAPEQ) 

COMPUTOR ALGORITHM TO CALCULATE LONGSHORE ENERGY FLUX FACTOR AND WAVE 
DIRECTION FOR TWO PRESSURE SENSOR ARRAY 

MAIN PROGRAM 


PRUGKAM I$ PRESENTLY BET UP TU TAKE A TIME SERIES OF 1024 POINTS IN MAIN 


OIMENSION C(Si2) 

DIMENSION FIR( 1024) oF IT(1024) oF2R(1024) Fel(ioed) 
DIMENSION SIGMA(512) oFMODSQ(S1R) pTHETACSI2) 
OIMENBSION C12(°S12)9012(512) 

OIMENSION w(ioe4) 

DIMENSION CG(S12) e8INTHB(S42) 

REAL MEAN] )MEANQ 

LOGICAL END 

DATA END/,FALSE,/ 

FURMAT($0(C2X0F 8,2) ) 


DEFINITIONS@FIXED VARIABLES 

KJEXPONENTIAL POWER DEFINING NUMBER OF TIME SERIES POINTS#(2%8K) 
S]SPACING BETWEEN WAVE GAGES (FEET) 

DELTTSTIME STEP BETWEEN POINTS IN AVERAGED TIME SERIES (SECONDS) 


BETA BANGLE DIFFERENCE BETWEEN WAVE GAGE ALIGNMENT AND SHORELINE (RADIANS) 


SLOPESSLOPE OF BEACH AT PUINT OF WAVE BREAKING 

GAMMAMSPECIFIC WEIGHT OF FLUID (LAS/FT#*#3) 

BSDISTANCE OF PRESSURE SENSORS ABOVE BOTTOM (FEET) 

GEACCELERATION OF GRAVITY (FEET/SEC##2) 

GBBRATIO BREAKING WAVE HEIGHT/DEPTH FOR LINEAR THEORY COMPUTATION 
OF WAVE HEIGHT 

GBPBRATIO BREAKING WAVE HEIGHT/DEPTH FOR LINEAR THEORY COMPUTATION OF 
WATER DEPTH GIVEN BREAKING WAVE WEIGHT 


DEFINITIONS©FLOATING VARIABLES 
AVGIBAVERAGE OF TIME SERIES) 
AVG2SAVERAGE OF TIME B8ERIES 2 
CCI) SWAVE CELERITY 
ClaCI)BCOSPECTRA OF SERIESIe2 
CBZBREAKING WAVE CELERITY 
CGC) SGROUP WAVE CE&LERITY 
CNTILCI)24096 POINT TIME SERIES BEFORE AVERAGING 
DEPTH@DEPTH OF WATER AT GAGE SITE FROM AVEHAGES OF GAGES 1! AND 2 
FLICI)BUNDEFINED/COMPLEX IMAGINARY PURTIUN OF TRANSFORM 
FIRCI)STIME BERIES DATA GAGEI/COMPLEX REAL PORTION OF TRANSFORM 
F21(1)BUNDEFINED/COMPLEX IMAGINARY PORTIUN OF TRANSFORM 
FeR(L)aTIME SERIES DATA GAGE2/COMPLEX REAL PORTION OF TRANSFURM 
FMUDSQ(T)e1IME SERIES AMPLITUDE MODULUS SQUARED 
HBBRREAKING WAVE HEIGHT 
TAC1)35000 POINT PATA GROUP AND TIME SERIES RECURD 
PLNEG3NEGATIVE CONTRIBUTIUN TO LONGSHURE ENERGY FLUX FACTOR 
PLNETENET LONGSHORE ENERGY FLUX FACTOR 
PLPOSzPOSITIVE CUNTRIBUTIUN TO LONGSHORE ENERGY FLUX FACTOR 
GQL@CT)BVYADSPECTRA OF SERIES 122 
RSSCALING FACTOR FOR SCALING UP ENERGY OF NONUSABLE 
PORTIONS OF DIRECTIONAL SPECTRA 
RATIUISRATIO OF ENERGY/WINDOWEUD ENERGY FUR GAGE 1 
RATTU2SRATIU OF ENERGY/WINDOWED ENERGY FUR GAGE 2 
REAL(I)#1024 POINT TIME SERIES AFTER AVERALING 
RN@RATIO OF GROUP WAVE CELERITY TO WAVE CELERITY 
RSHFRUSPERCENT OF ENERGY BEYUND SPACIAL ALJASING FREQUENCY 
RSLFRQRPERCENYT OF ENERGY BELUW LUW FREQUENCY CUTOFF 
RSUODSPERCENT OF INCUHERENT ENERGY : 
SHFREQ@SSUM UF ENERGY WITH FREQUENCIES ABOVE SPACJAL 
ALIASING FREQUENCY CUTOFF 
SHG2BSUMMATION OF ONSHORE ENERGY FLUX 
SIGMACI)BRADIAL FREQUENCY i 
SLFREQ@SUM OF ENERGY WITH FREQUENCIES BELOW LOW FREQUENCY CUTOFF 
SUDD@SUM OF ENERGY WITH FREQUENCIES HAVING INCOHERENT WAVE DIRECTION 
SUMEN@SUM UF ENERGY 
SUMISSUM OF SQUARES OF TIME SERIES 1 WITHOUT AVERAGE 
SUM2MSUM OF SQUARES UF TIME SERIES 2 WITHOUT AVERAGE 
SUMF{e8IIM OF SQUARES OF TIME SERIES 4 WITH AVERAGE 
SUMF228UM OF SUUARES UF TIME SERIES 2 WITH AVERAGE 
TEWAVE PERIOD ES 
THETACI)@WAVE DIRECTION IN RADIANS 


Figure 2. Listing of main program. 
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73 


80 


6s 


90 


93 


120 


125 


135 


aan 


aonanarann 


oan 


aonn 


29 


30 


TMETABSBREAKING WAVE ANGLE 

WSUM{@SUM OF SQUARES OF DATA WINDOW MODIFIED TIME SERIES 1 
WSUM2@8SUM OF SQUARES OF DATA WINDOW MUDIFIED TIME SERIES 2 
Pla3,141$9265 

TWOP La2.O*Pl 

Keio 

Nae, 88K 

$260,0 

DELTTs{.00 

BETAS1,$708 

SLUPE&O,U§ 

GAMMAz64,0 

bs2.5 

MaNe]l 

Nb2sN/2 

GB20.78 

6232.2 

GRPsU.78 


HIGH FREQ CUTOFFAa3,0 SEC 

SPACIAL ALIASING CUTOFF8304 SEC 

NLOWSLOW FREQUENCY CUTOFF NUMBER 

NYFR@HIGH FREQUENCY CUTOFF NUMBER 

NSALFRSSPACIAL ALIABSING FREQUENCY CUTOFF NUMBER 
FREQUENCY CUTOFF NUMBER®TIME SERIES LENGTH/CUTOFF PERIOD 


NLOWS50 
NYFRa3aa 
NSALFRa301 
NSMIBNSALFRe1 
CONTINUE 


INITIALIZING VALUES 
SLFREQGs0,0 

8UDD#80,0 

SHFREOQ#0,0 
SUMENEO 0 

SUMi#0, 

SUM2z0, 

SUMFigy, 


SUMFez0. 
WSUM150,0 
wSUM280,0 
AVG180,0 
AvG280,0 


PLPOS80.0 
PLNEGa0.0 
PLNET@#0.0 
8HG280,0 

DUO 29 ImieN 
FLl(1)#0.0 
F21(1)20.0 
CONTINUE 

00 30 Tal»ND2 
FMUDSG(1)80.0 
CONTINUE 


THIS PORTION OF PROGRAM READS IN WAVE PRESSURE VALUES INTO FIR,F2R ARRAYS 
AND ASS8URES MATCHING DATE GROUPS FUR DIRECTIUNAL WAVE ANALYSIS OF TWO 
G 
CALL BUF ENGAGE, sMONTHL pMDAYIeMTIMEZ9FIR 9 IDATES 9 END) 
IFCEND) GU TO 1 
CALL BUF(MGAGE2sMONTH2eMOAYBeMIIME2,F2R eIDATE2sEND) 
IF(END) GO TO 1 
IF(IDATE1,€0,%DaT22) GO TO 120 


Figure 2. Listing of main program.——-Continued 
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140 


145 


190 


185 


160 


165 


175 


180 


185 


190 


195 


200 


aan 


ann 


onn 


(ol eit ei ed 2) 


ue 


4} 


Ag 


43 


97 


BACKSPACE 9 

GU TU 110 

CONTINUE 

IF(MGAGE1,NE,311) CALL SWITCH(MGAGEY »MGAGE2@oF{R oF aR ) 
WRITE (60426) 

FURMAT(//0( GAUGE NO, looXe (MUNTH (9 9X9 (DAY Le BXe (TIME () 
WRITE COo11)MGAGE1LoMONTHS oMDAYIOMTIMES 

WRITE (6011) MGAGE2eMONTH2 eMDAV2OMTIME? 
FURMAT(I793(SxX0977)) 


THIS PORTION OF PROGRAM CALCULATES WATER DEPTH AT WAVE GAGES AS WELL as 
AVERAGES AND SUM OF SQUARES UF TIME SERIES 

DO 4@ Jmi_N 

AVGIBAVGI+FIR(T) 

AVG2BAVG24F 2R(1) 

AVGIBAVGI/FLUAT(N) 

AVG2BAVG2/FLOAT(N) 

DEPTNa( AVG) +AVG2) /20%B 

CALL HFC( DEPTH, So DELTT)»NoNSALFR) 

OU 41 Taton 

FIRCL)@FIR(I)eAVGI 

FeR( I) sFaR(1)eAVG2 

SUMIZSUMIFFIR(T)*%2, 

SUM2ZSUMO+FAR( 1) **%2, 

CONT INUF 

SUMLESUM] /FLUAT(N) 

SUMABSUMA/FLOAT(N) 


THIS PORTION OF PRUGRAM APPLIES DATA WINDOW TO TIME SERTES@=@DATA WINDOW 
VALUES ARF REPRESENTED BY WC1) 

DO B9 Jalen 

WOT) 90.5% C1 0SCUSCTWUPIMFLUATCIY/FLUAT(N))) 


FIRCLBCFIRCI) yew(1) 
FaRCI)a(FaRCT) powcr) 
CUNTINUE 


THIS PORTIUN OF PROGKAM COMPUTES SUM OF SQUARES OF DATA WINDOW MOOIFIED 
TIME SERIES AS WELL AS RATIO OF PRE WINDUWED ENERGY TO WINDOWED ENERGY 

DU 43 Jaton 

NSUMLaWQUMIFFIR( 1) *#2, 

WSUM@BWSUMO+FAR(IT)**2, 

CONTINUE 

WSUM}BWSUM] SPLOAT(N) 

WSUM2EWSUMA/FLOAT(N) 

RATIUIaS9UM{ /w8UMI 

RATIU2sSUM2/WBUMe 

CALL FETC FYRoFiToKo0) 

CALL FFT(FeRoFeteKe0) 

MEANJ@FIR() 

MEAN@EFaR(1) 


THIS PORTION OF PROGRAM CALCULATES CO AND YUAD SPECTRA VALUES, AS WELL A8 
WAVE ANGLE 7O SHORELINE AND ENERGY CONTRIBUTIONS OF EACH FREQUENCY, 
BREAKING WAVE HEIGHT AND BREAKING WAVE CELERITY ARE ALSO 
CALCULATED IN THIS SECTION 

Ta1 

DO 97 JB2oN 

FIR(I)sFIR( J) 

FLLCL)BFIT(J) 

FeR(1)sFeR(J) 

F2l(1)eFal(J) 

Tele) 

CONTINUE 

DO 96 yalom 

SUMFIBSUMFYSFIR( I) #*R oF LIC 1) 9%2, 


Figure 2. Listing of main program.——Continued 


16 


205 SUMF2aSUMFa+FaR(1)¥#2,¢F2I(1)#¥2, 
90 CUNTINUE 
SUMF 1SSUMF{*MEANIF#2, 
SUMF 2BSUMF 2+MEANG&#2, 
WRITE (60289) 

210 289 FURMAT(//e7Xe lI felOXe (GSIGMACT) Lol ixe (FMDBQ(T) () 
DU 99 [TmieND2 
Ci2(l)yeFiR(I)#FaRClyeFil(l)sFei(l) 
GQ12C1)aFiR( I) *Fal(LyeFeR(I)#Fil(l) 
SIGMAC]) SFLOATC IT) #THUPT/C(CFLOAT(N) #DELTT) 

ais TSTWUPT/8IGMACT) 

CALL wVLEN(DEPTHeTeXKH) 
XKBXKH/DERTH 

IF (C12(1) LE.0.00000000!1) GO TU 98 
POE(1e/(KKHS)) RATANCULECT) /C12(1)) 

220 IF(ABB(PD),GT,1.0) GU TU 93 
THETA(I)s=@ACOS(PD)¢BETA 
GO TU 92 

95 THETACI)80,0 
GU TU ga 

223 93 THETA(I)20,00001 
FMODSOCT)BFIR(I)**2,tF11(1)**20 
XKBEXKH*B/DEPTH 
XKP=CUSH(XKB) /COSH(XKH) 
FPMUDSO(1)sFMODSQO(1) /(XKP8#2,) 

230 FMODSQ(I)sFMODSQ(T)¢#RATIOL 
SUDDSSODD+FMUDSQA(1) 

WRITE( 95105) Tp SIGMACI) oFMUDSUCI) 
105 FORMAT(3Xo TS pSxXoFigeoe7XoF 1 Reb) 
92 CUNTINUE 

235 FMODSQ(I)aFIR(I)##aetFi1(1) 20 
XK bSAKHSB/DEPTH 
XKPSCOSH(XKB) /COSH(XKH) 
FMODSQ(T)SFMODSQ(1) /(XKP#%2,) 
FMUDSQ( IT) sFMODSQ(I)*RATION 

240 RN®0 e580 1 +2 e*XKH/SINH(2,*XKH)) 
CCCI) ag IGMA(l) *DEPTHERN/XKH 
CCI) SCGCI)/RN 
HG2e(CG(1)#2,08FMUDSU( 1) *COS(THETA(S))) 

Cc SHG2SONSHORE ENERGY FLUR 

2a5 : ShG288HG2+HGe 
SUMEN@GUMEN+F MODS U( I) 
IFCI.GE,NSALFR) GO TO 79 
GO TO ¥B8 

79 SHFREQsSHFREQs+FMODGQ(T) 

250 78 CONTINUE 
TF(I,LEeNLOW) GO TO 77 
GO TU 76 

77 SLFREG=8LFREQ+FMUD8Q(T) 
To CUNTINUE 
255 1F(1.GE,NYFR) GO TO 999 
99 CONTINUE 
999 CONTINUE 
SHG2ESHG2e2, 
WRITE (60351) 

260 351 FORMAT (//oOXo (1 tolaXe (SIGMACI) besiXe (PCT (oloXe (THETACS) Cl) 
DU 4B Jai yNDe 
IFCI,GE.NSALFR) GO TO gu4 
PCTSFMOD§0(1) /SUMEN 
IF(PCT,GE.0.0a5) GO 10 49g 

265 GO TU 4@ 

49 WRITE (6050) ly 8IGMACI) oPCTeTHETACT) 
SO FURMAT( 3X9 1S93(3X0F160,.8)) 

4B CUNTINUE 

44 CONTINUE 


Figure 2. Listing of main program.—-Continued 
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270 ; HBm(B,#8,4)*(8HG2*%, 4) e(GH/G) 8,2 
CBa(G*HB/GBP) ae 5 
57 CONTINUE 


THIS PORTION OF PROGRAM MUDIFIES WAVE GAGE ANGLES TO BREAKING WAVE ANGLES 
AND COMPUTES LONGSHORE ENERGY FLUX FACTORS 

DU 91 Ieienb2 

TF(I.GE.NSALFR) GO TO 996 

SINTHBCIT)SS8IN(THETACL))*CB/C(I) 

THE TABBASINC(SINTHB(I)) 

280 KKRSB((LoeSINCTHETACT))®#20)/( 1 o@SINTHB(1)*¥2,))**,5 
XKSS§CG(1)/CB 
FMUDSO(1) BFMODSO(1) *XKRS#XKSS 
IFC THETA(1) oLE.020) bU TU 87 
PLPUSEPLPOS+GAMMA¥SIN(2,*THETAB )*#CB*FMUDSQ(T) 

285 GO Tu Og 


B7 PLNEGEPLNEG+GAMMA*SIN( 2, *THETAB 2*®CB*FMODSO( 7) 
86 CUNTINUE 
PLEGAMMASCBSSIN(2,*THETAB )*FMODSQOC(T) 
PLNETSPLNET¢PL 
290 IF(I,GE.NYFR) GO TO 998 
91 CONTINUE 
998 CUNTINUE 
RSODUSSODD/SUMEN 
KSHFROSSHFREO/SUMEN 
295 RSLFRUSSLFREQ/SUMEN 
RTUTSRBUDD+RSHERO 
REL/( 1 ,eRTOT) 
PLPUSBPLPUS#R 
PLNEGSPLNEG#R 
300 PLNETSPLNET#R 
WRITE C0201) NSALFR 
201 FORMAT(//o[ NSALFRwle2uXelo) 
WRITE (60200) DEPTH 
200 FORMAT(( DEPTH UF WATER AT GAUGE SITES (oF10,1) 
305 WRITE (60100) AVG1 9 AVGa 
100 FORMAT( ( AVG1a(oF14.399X0 [LAVG2S (oF 44,3) 
WRITE (60170) SUM{ 9 SuM2 
170 FURMAT(( BSUM{e (oF 11, 309Xe (8SUMAE (oF 44,3) 
WRITE (60111) WSUML 9 WSUMO 
310 111 FURMAT( [ WSUM{ a (oF 100399X0 (WSUM2B (oF 10.3) 
WRITE (COoL12)RATIVUL,»RATIO’ 
112 FORMAT( ( RATIOS (oF 9039 9X0 (RAT LOZ@ (9 F953) 
WRITE (6039) SUMEN 
39 FORMAT( ( SUMENS [yo PXoF 13.8) 
R15 WRITE (CO0104) HB 
104 FOKMAT( ( BREAKING WAVE HEIGHT HBe8()6xX0F10,2) 
WRITE (60108)C8 
108 FURMAT( { BREAKING WAVE CELERITY Cam (o4xoFl0,2) 
WRITE( 60106) RSODD»RSHFROVRSLFRO 
320 106 FORMAT(( R§O0DDM (oF 11,49 8X0 IRSHFROS (oF 10,49 8X9 (RSLFROB (9 F 10,4) 
WRITE (694$03)PLPOS0P LNEG 
103 FURMAT([( PLPOSB [oF il o4,8Xo (PLNEGw ly fll ,4) 
WRITE (60109) PLNET 
109 FORMAT(( PLNETa (oF 11.4) 
35 GU TU 410 
1 CONTINUE 
8TUP 
END 


canon 


275 


Figure 2. Listing of main program.--Continued 
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(1) Input data for this program are in the form of digital 
magnetic-tape records of 4,100 values. The first 4 values of the 
records are the gage number, month, day, and time of the observa- 
tions; the last 4,096 values are the time-series pressure values of 
the wave gage. In the present program the wave gage pressures are 
stored in thousandths of a foot (head) water at 0.25-second inter- 
vals. Subroutine BUF reads time-series data into array CNITL, 
where it is averaged to provide 1,024 time-series values of At = 1 
second spacing. Units are also divided by 1,000 to convert values to 
feet (head) of water. 


(2) The date groups of record 1 and record 2 are compared to 
ensure that times of records are simultaneous; if the times are not, 
the program searches the record file until this condition is met. 
The two records are than checked for proper sequence to ensure that 
gage 1 is analyzed first. Subroutine SWITCH switches arrays if they 
are not in proper order. 


(3) Each of the two 1,024 value time series is then analyzed for 
average values which are printed out along with the average depth of 
water at each gage site. The average value of each of the time- 
series records is again averaged and is added to the height of the 
gages above the bottom to obtain the water depth: 


AVERAGE 1 + AVERAGE 2 | 


EPTH = 
DEPT 5} 


B 


in which AVERAGE 1 is the average of time series 1 = a,(0), AVERAGE 2 
the average of time series 2 = a, (0), and B the height of sensors 
above the bottom. 


An option to apply a weighting function w(j) (= W(1) in program) 
has been incorporated before the FFT subroutine is called. Im this 
particular program a cosine bell weighting function has been incorpo- 
rated. If the data window option is selected, the two time-series 
data records, which are read into FIR and F2R arrays, are multi- 
plied by the following weighting function (cosine bell) 


w(j) = aE - cos (2)| 


where j is the time step number and N the number of data points 
in series. If no weighting function is desired in analysis set 
w(j) = 1.0, which is the “box car" weighting function. 


As the cosine bell function reduces the total energy content of 
the waves, the final energy obtained from the FFT must be rescaled up 
to the proper value. This is accomplished by scaling up the time- 
series pressure values by the ratio 


Unwindowed energy <p2> 
R = = 
Windowed energy <p'2> 


as discussed in equation (4). 
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(4) Cospectra and quad-spectra of the gages are computed using 
the following relationships (note in computer program index, I is 
used for frequency counter, n): 


Cospectra = C12(1) = FIR(1)*F2R(1L) + FII(1)*F21(1) 


Quad-spectra = Q12(1) = FIR(1)*F2I(1) - F2R(1)*F1I1(1) 


in which FIR and FII are the real and imaginary parts of complex 
transforms of time series 1: F2R and F2I are the real and imagi- 
nary parts of complex transforms of time series 2. 


(5) Wave angle is calculated in accordance with equation (19). 


peers feel iy 1 alan) 
Oa) SO & arcosine Een cope eae Heal 


where k(n) is the wave number calculated via linear wave theory, 2 
the spacing of gages, and £8 the difference in alinement of gages 
and shoreline in Figure 1. 


Due to energy leakage problems in spectra, impossible wave angles 
can result [wave angles with (1/k(n)& arctan Q12(n)/C12(n)) greater 
than 1.0]. When this happens, energy is lumped into a separate 
category for later scaling up of the longshore energy flux. 


(6) The high frequency cutoff in this particular program has been 
set at 2.09 radians per second, which corresponds to a period of 3 
seconds or NYFR = 342. This value can be reset in the main program 
by adjustment of NYFR where 


and N is the number of data points in time series, At the spacing 
in time of data points, and Tur the high frequency cutoff period. 
The spatial aliasing frequency is computed in subroutine HFC. 


Energy between the spatial aliasing frequency and the high fre- 
quency cutoff is put into a special category and used to scale up the 
final longshore energy flux. 


(7) Each frequency contribution to the onshore energy flux is 
calculated for the gage site location as follows: 


Onshore energy flux = 2y EXGyili C,(n) cos [0(n) ] 
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|F,, (2) | = modulus of the complex amplitude spectra of wave 
elevation above mean surface at gage site 

Cy (n) = group wave speed at gage site 

O(n) = © = angle of wave direction (see Fig. 1) 

Y = specific weight of seawater 


The onshore energy flux is then summed to obtain the total onshore 
energy flux. In the program, onshore energy flux/y = HG2. 


(8) Breaking wave height at the shoreline is determined from the 
mean square onshore energy flux via a linear theory wave transforma- 
tion formula which can be simplified to 


N/2 cB\0:2 
LS as 16|F (nm) |? Cc, (n) cos O(n) o«( =) 


where GB is the wave height-to-water depth ratio at breaking and g 
the acceleration of gravity. 


The choice of GB is up to the user although a value of GB = 
1.42 has been found by Komar and Gaughan (1972) to best fit wave tank 
data for breaking wave heights for monochromatic waves. In the pres-— 
ent program, GB has been set equal to 0.78 but can be readily 
changed. 


The breaking wave water depth is calculated from the equation 


Bb 
dy 


= GBP 
where d is the wave breaking water depth and GBP the ratio of 
wave height to water depth at breaking. 

In this case a different value of the ratio of breaking wave 
height to water depth can be used in the program for obtaining the 
proper water depth. An assumed value of GBP = 0.78 from the solitary 


wave theory in the SPM is used. 


Linear wave celerity is assumed and breaking wave celerity is 


estimated as 
H. \0°5 
; b 
a, « (c%) 
2 (. 


The breaking wave height and celerity calculated in this approach 
apply to all frequencies. 
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(9) Frequency—by-frequency modification of wave angles is made 
assuming linear wave theory, Snell's law, and parallel bottom con- 
tours offshore. The breaking wave angle, O,(n), is calculated from 


C,(n) sin 0,(n) 


C 


0,(n) = arcsin | 
r 


where the subscript r refers to the reference gage location. 


(10) Longshore energy flux is calculated for each frequency 
component (except the special cases discussed in Sec. II) using the 
equation 


Peg(n) = y|F,(n)|* C.4(m) sin 20,(n) 


Cg 
and is summed up to obtain a net longshore energy flux. 


(11) The value of the net longshore energy flux is multiplied by 
a factor R which scales up the total energy in the spectrum (below 
the high frequency cutoff). The equation for scaling factor R is 


1 


R= (7 — RTor) 


where RTOT = RSODD + RSHFRQ when RSODD is the percent of energy in 
low frequency bands for which impossible values of the cosine func-— 
tion are calculated, and RSHFRQ is the percent of energy between 
spacial aliasing frequency and high frequency cutoff. 


The final result of analysis of the two gage records for the 
net longshore energy flux PLNET is printed out, as well as specific 
frequencies for which impossible directional results occur and fre- 
quencies at which more than 2.5 percent of the total wave energy is 
found. 


IV. SUBROUTINE DOCUMENTATION 


FFT Subroutine. 


The sampled time function, f(j), will be expressed as 


N-1 
£(j) = ) F(n) exp(inw, jAt) 
n=0 
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in which 


1 record length aE NAt 


t. = jAt = a discrete time where j is the integer time step 


F(n) = a(n) —- ib(n) 


N = N N 
pea ail ad) it omad? 1 


Ne o\e N N 
b (55a) = Bl Gea) eee; 
a(0) = mean of sampled record 


b(F)= 0 


Because negative indexes are not readily handled by most FORTRAN compilers, 
the summation extends over the interval 0< n< N- 1, rather than over the 
symmetric interval - N/2 $n $ N/2. From the definition of the coefficients 
above, it is clear that the coefficients a(n) and b(n) for n > N/2 contain 
no additional information. 


b(0) 


The inverse relationship completing the FFT pair is 


N 
1 
F(n) =— ) £(j) exp(-inw, jt) 
N; " 
j=1 
As an enumeration of the complex FFT coefficients, suppose the series of 8 
values is considered, N = 8. The coefficients would be 


F(O) = a(0) 
F(1) = a(1) - ib(1), F(7) = a(7) - ib(7) = a(1) + ib(1) 
F(2) = a(2) - ib(2), F(6) = a(6) - ib(6) = a(2) + ib(2) 
F(3) = a(3) = ib(3), F(5) = a(5) — ib(5) = a(3) + ib(3) 
F(4) = a(4) 


This pattern prevails for all sets of FFT coefficients, regardless of the 
value of N. Both F(0) and F(N/2) are real and, as noted previously, the 
coefficients F(n) for n> N/2 really contain no additional information. 
The FFT subroutine used here requires that the number of data points, N, 
provided be an integral power of 2, i-e., 


n= 2K 


where K is an integer. Thus analyses of 512, 1,024, or 2,048 data points 
(K = 9, 10, 11) would be suitable with this subroutine. 
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The following two requirements are satisfied in the FFT subroutine. 


(a) By operating on the sampled function, obtaining the F(n) 
coefficients and carrying out the inverse FFT (FFT !), the original 
time function is recovered. Schematically, 


(b) The mean square of the sampled time function is equal to the 
sum of the squares of the moduli of the FFT coefficients, F(n), i.e., 


l N N-1 
=) f£Gle =) Seq) 
Nga n=0 


a. Calling Statement: SUBROUTINE FFT (FR, FI, K, ICO) (see Fig. 3). FR, 
FI = real and imaginary coefficients in 


F(n) = FR(n) —- iFI(n) 
K = power of two (i.e., N = ay 
ICO = control whether FFT or (FFT)~! 


operation is desired if 


= 0 > FFT 


Ico E 
1 > (FFT) 1 


When entering the subroutine, FR is the time sequence f(j) and FI is 
arbitrary. When exiting the subroutine, FR and FI are the real and imagi- 
nary parts of the complex transform, respectively; e.g., input is 


K=5 
ICO = 0 
Lisi 2n(jAt) 4m(j At) 
£(j) = 1.0 + 2.0 cos “aya 3.0 cos “Tn aps 
o a uC) .. 4n(jAt) 
0.6 sin aay 1.4 sin gia: 
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Db. 


10 


15 


20 


25 


30 


35 


40 


a3 


50 


Figure 3. 


Data Input to Program. 


FAST FOURIER TRANSFORM SUBROUTINE 
SUBRUUTINE FFIT(FReFIe Ky ICO) 
DIMENSION FR(1)9FIC1) 
NSQREK 
IF (ICO,&G,0) GO TO 10 
DU B Twin 

6 FIC T)SeeF ICT) 

10 CONTINUE 
MREO 
NNENe]} 
DUO 2 MaiegNNn 
Lan 

1 L3l/2a 
IF(MR+L.GT,NN) GO TO } 
MREMUD(MReL) +L 
IFCMR, LEM) GO TO 2 
TREFR (M41) 
FROM¢1) BFR(MRO1) 
FR(MR¢{)87R 
TISFICM41) 
FI(M@1)eFI(MRe1) 
FICMRe{)aTy 

2 CONTINUE 
Lal 

3 IF(L.GE.N) GO TO 7 
ISTEPwaeL 
ELSL 
DO 4 Msiol 
ABS.14159ROS3S5#FLOAT(1=@M) JEL 
wRaCUS(A) 
wIBSIN(A) 
OO WY ImMyNe ISTEP 
valet 
IFCICU,EG.1) GO TO 11 
TREWR*ER( J) eoHTOFI(J) 
TISWReE T(J) +H eFR OJ) 
GO TU 42 

11 TREWR#FR( J) OWT SFI CJ) 
TISWREFI(J)oHTePR( J) 

12 FR(JJSFR(I)°TR 
FICJ)SFICL)@TI 
FRCI)8FR(I)¢TR 

G FICI)@FICI)¢TY 
L=IsTeP 
GO TU 3 

7 CONTINUE 
ANEN 
IF(ICO,€0,41) GO TO 6 
DO § Ie19Nn 
FRCI)SFRCI)/AN 

S FICI)seFICI)/AN 

6 RETURN 
END 


_o0 


Listing of FFT subroutine. 


£(j) values at 6-000 5-080 3.750 
At = 1 second 0.590 -0.829 -1.900 
intervals -2.600 —2.215 -1.451 
(32 values) 0.562 1.445 2.034 
2.000 1.391 0.513 

-1.390 —2.054 =D e322 

-1.400 =—0.257 1.188 

4.238 5.438 6.189 

FR = 6-000 5.080 3.750 
(32 values) 0.590 -0.829 -1.900 
—2.600 I) 72 1N5) =1].451 
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2-184 
—2.506 
-0.465 

2.229 
-0.475 
-2.109 

2.755 

6-386 

2.184 


-0.465 


0.562 1.445 2.034 2.229 


2.000 1.391 0.513 -0.475 

-1.390 -2.054 —2.322 -—2.109 

-1.400 -0.257 1.188 2.755 

4.238 5-438 6.189 6.386 

FI = 0.000 0.000 0.000 0.000 
(32 values) 0.000 0.000 0.000 0.000 
0.000 0.000 0.000 0.000 

0.000 0.000 0.000 0.000 

0.000 0.000 0.000 0.000 

0.000 0.000 0.000 0.000 

0.000 0.000 0.000 0.000 

0.000 0.000 0.000 0.000 

c. Calling Statement: FFT (XR, XI, 5, 0). 

Output: 1.000 1.000 1.500 0.000 
a(n) coefficients 0.000 0.000 0.000 -0.000 
(32 values) -0.000 -0.000 -0.000 -0.000 
-0.000 -0.000 -0.000 -0.000 

-0.000 -0.000 -0.000 -0.000 

—0.000 -0.000 -0.000 -0.000 

-0.000 0.000 0.000 0.000 

0.000 0.000 1.500 1.000 

b(n) coefficients 0.000 -0.300 -0./00 -0.000 
(32 values) -0.000 -0.000 -0.000 -0.000 
-0.000 -0.000 -0.000 -0.000 

-—0.000 -0.000 -0.000 0.000 

0.000 0.000 0.000 0.000 

0.000 0.000 0.000 0.000 

0.000 0.000 0.000 0.000 

0.000 0.000 0.700 0.300 


At (time step) = 1 second in above example. 


2. HFC Subroutine. 


This subroutine resets the spatial aliasing frequency cutoff to a higher 
frequency than would be the case for normal incidence of waves to gage pair. 
In the present version of this subroutine, it has been assumed that the maxi- 
mum angle which the wave crests can make with the gage pair axis is 45°. The 
spatial aliasing criteria are expressed in Figure 1, where for proper resolu- 
tion of wave direction the following criteria must be met 


Rcos [Q(n) - 8B] <3 


k(n) £cos [O(n) - 8] < k(n) 


NS} 
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The proper spatial aliasing frequency to correspond with the spacial aliasing 
wave number cutoff is found from the normal wave dispersion relationship. 


Calling Statement: HFC (DEPTH, S, DELTT, N, NSALFR) (see Fig. 4). 


DEPTH 


DELTT 


NSALFR 


10 


depth of water at gage site (from main program) 


spacing of wave gage pair (from main program) (= & in text) 


time-step increment between values in time series analyzed 
(from main program) 


exponent of 2 describing number of time series values 
(from main program) 


integer number for spatial aliasing frequency cutoff 


om oaaon 


SUBRUUTINE HFC(HIGH FREQUENCY CUTOFF/SPACIAL ALIASING FREQUENCY) 
RESETS ALIASING CUTOFF TO HIGHER FREQUENCY 
HASED UN ASSUMED MAXIMUM WAVE ANGLE 

SUBRUUTINE HFC(DEPTHsS»)DELTToN¢NSALFR) 

SPACIAL ALIASING ASSUMES WAVE ANGLES LESS THAN 45 DEGREES 

XK829,14159/0.707 

KKH@AK§aDEPTH/S 

§1GSUa32,2¢(XKH/DEPTH) #TANH(XKH) 

SIGHF aBORT(S81GSQ) 

RECLNSFLOAT(N) SDELTT 

NSALFRaSIGHFSRECLN/60283 

RETURN 

END 


Figure 4. Listing of HFC subroutine. 


3. SWITCH Subroutine. 


This subroutine is set up to interchange time-series data arrays in the 
instance when gage 2 data are processed before gage 1 data (see Fig. 5). If 
the first gage record processed is not equal to the appropriate number of the 


gage, 


as specified in main program, data arrays of first and second gage 


records are interchanged. 


10 


15 


4. WVLEN Subroutine. 


aan 


SUBRUUTINE SWITCH EXCHANGES LOCATIONS OF TIME SERIES DATA TO ASSURE 


GAGE1 I8 STORED IN FIRST ARRAY AND GAGE2 IN SECOND 


SUBRUUTINE SWITCH(M1oMaeFIReFeR) 
DIMENSION FIR( 1024) oF2R( 1024) oF 3R(1024) 
M3aM{ 

Misha 

MaeMy 

DO 10 Ieiy102u 

F3R(I)@FiR(I) 

FIR(CI)erFaRr(l) 

FeR( 1) eF3R(I) 

CONTINUE 

RETURN 

END 


Figure 5. Listing of SWITCH subroutine. 


This subroutine accepts wave period and water depth as input and calcu- 
lates wave number as output via a Newton-Raphson iteration. 
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Calling Statement: WVLEN (DPT, PER, XKH) (see Fig. 6). 


DPT = water depth (from main program) 


PER 
XKH 


wave period (from main program) 


wave number * water depth (calculated in subroutine) 


WAVE LENGTH ITERATION SUBROUTINEee THIS SUBROUTINE CALCULATES WAVELENGTH 
VIA NEWTON@RAPHSON ITERATION USING PERIOO,WATER DEPTH INPUT 
PEKBWAVE PERIOD 
OPTSWATER DEPTH 
XKHOWAVE NUMBER®WATER DEPTH 
BUBRUUTINE WVLENCOPT 9 PER» XKH) 
KKHOS(6,2631853/PER) #8 2¥DPT/3202 
IF (XKHOS6,3)2e1ok 
10 1 XKABXKHO 
GO 70 9 
XKHBSURT(XKHO) 
SHESINH(XKH) 
ChHECUSH(XKH) 
15 EP MBX KHUCKKHSSH/CH 
SLOPEmeXKH/CH#*2—8H/CH 
DXKH@eEPS/ SLOPE 
IF (ABS (DXKH/XKH)00,0001)99994 
4 XKHEMKHeDKKH 
20 GU TO 3 
9 CONTINUE 
RETURN 
END 


fetete lel ata] 


i 


Figure 6. Listing of WVLEN subroutine. 


5. BUF Subroutine. 


This subroutine is set up to read in wave gage files from magnetic tape. 
The data records consist of arrays of 4,100 values, the first four of which 
are the gage number, month, day, and time of wave record. The remaining 4,096 
values represent pressures in thousandths of a foot (head) water. The data 
are returned to main program aS a wave gage number-date series and a time 
series of 4,096 values of pressure in feet (head) of water. Two records are 
processed in one pass. 


Calling Statement: BUF (MGAGE, MONTH, MDAY, MTIME, CNTL, IDATE, END) (see 
leo 7) 0 


MGAGE = number of gage (read from tape) 


MONTH = month of observation (read from tape) 

MDAY = day 

MTIME = time 

CNTL = control array of 4,096 pressure values in feet (head) of water 
returned to main program 

IDATE = summed time group for time comparison between gages 

END = logical end 
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SUBRUUTINE BUF READS IN WAVE GAGE DATE INFU AND TIME SERIES DYNAMIC 
PRESSURE VALUES IN FEET HEAD UF WATER 
THIS SURROUTINE KEADS 4096 TIME SERIES VALUES AND AVERAGES TO OBTAIN 
1024 VALUES FOR MAIN PROGRAM ANALYSIS 
MGAGESGAGF NUMBER 
MONTHSMUNTH OF RECORDING 
MDAYZnAY OF RECORDING 
MTIMESTIME OF RECORDING 
REALM@ARRAY OF AVERAGED TIME SERIES VALUES 
SUBRUUTINE BUF (MGAGEsMONTHeMDAY oMTIMESREAL? IDATE END) 
DIMENSION CNTL(4096)91A(5000) 
DIMENSION REAL(1024) 
LOGICAL END 
15 DO 12 J™104096 
CNTL(J)80.0 
12 CONTINUE 
BUPFER IN(9e1)(TA(1)01A(9000)) 
IFCUNITC9) 10920030 
20 20 PRINT iio (IACI) o131049) 
11 FURMAT(( PARITY ERROD ON(e4I7) 
10 CUNTINUE 
MGAGEBIA(1) 
MONTH@JA(2) 
25 MDAY@IA(3) 
MTIMESTA(4) 
IDATEBTAC2) +I ACI) PIAC4) 
Key 
DO 25 Jmie4096 
30 KBKe] 
CNILC J) BIACK) 
25 CNTLCJ)SCNTL(J)/10000 
00 26 J84088.4096 
26 CNIL(J)@CNTL (4087) 
335 Js} 
DO 27 Leieoi024 
REAL(L)@(CNTLOEJ)OCNTL CUO) SCNTEC Joa) eCNThL(ds3))/40 
JaJeq 
27 CONTINUE 
40 RETURN 
30 END#,TRUE, 
RETURN 
END 


(elelalelel of of of ei a) 
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Figure 7. Listing of BUF subroutine. 
V. SAMPLE OUTPUT 


Three examples of output are presented for different dates for the wave 
gage pair at Channel Islands Harbor (Fig. 8). The year the data was taken was 
1975. 


The first set of frequencies lists amplitude modules squared of wave data 
having impossible direction results. The sum total of this energy (in decimal 
percent) is listed as the quantity RSODD in the variable output at the bottom 
of the output. Im the case of the wave data taken on 7-26-1600, the inco- 
herent data amounted to 0.004 (0.4 percent) of the total energy in the wave 
record. 


The second set of frequencies listed provides the wave direction for the 
frequency bands having a significant part of the energy (22.5 percent). In 
the case of the wave record taken on 7-26-1600, it is seen that the wave angle 
is reasonably consistent from the frequency-to-frequency band and is approxi- 
mately 0.70 radian (40.1°). 


The variable list provided at the bottom of the sampled output gives 
values of most importance in the analysis of wave information for longshore 
energy flux. The longshore energy flux output is in pounds per second units; 
the output in the first example is 89.23 pounds per second. 
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Example 1 


GAUGE NO, MONTH DAY TIME 
341 y 26 1600 
3ia 7 26 1600 
T SIGMa(1) FMpsa(1) 
3 0016408 20000288 
4 0084544 0000014 
5 0030680 2000087 
7 2042954 0000036 
9 0055223 0000012 
10 0061359 2000081 
41 0067495 0000048 
14 0085903 0000005 
24 0147262 0000187 
25 0153398 0000041 
2? 0165670 0000002 
29 0177942 2000099 
30 0184078 : 0000040 
32 0196350 2000066 
33 0202485 2000029 
4a 0257709 2000014 
us e27O4i7 2000041 
65 0398835 2000093 
I SIGMACT) PCT THETA(T) 
o? 041110685 003904604 ©70276903 
68 e413 72402a77 004032988 077998085 
73 044792439 00P8312S1 068437028 
74 04540583) 010395194 069738613 
75 046019424 006890760 069498090 
78 04786020) 002500035 o7191 2354 
79 048473793 204487282 06379139 
NSALFRe 201 
DEPTH OF WATER AT GAUGE SITES 23.2 
AVG\s @yo4il AVG28 19,999 
8UM;s 0229 SUM2s 2234 
WSUMis 2084 WSUM@e 0086 
RATIO{e 2e730 RATIO2s 2.729 
SUMENSs 018433938 
BREAKING WAyE HEIGHT Hes 3.03 
BREAKING WAVE CELERITY CBB 11.18 
R&8Q00Da 00040 RSHFROs 02U13 RSLFRus 00170 
PLPOSs 9y 64a) PLNEGe 25,4150 


PLNETa 89.2271 


Figure 8. Three examples of output for wave gage pair at Channel 
Islands Harbor. 
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Example 2 


GAUGE NO, MONTH DAY 
311 T 20 
312 7 26 

I SIGMA(1) 
i 0006136 

2 v0s2e72 

a} 0038408 

b) 2030680 

8 0049087 

9 2055223 
10 0061359 
11 00607495 
13 0079767 
14 2085903 
15 2092039 
16 0098175 
18 elsouu7 
19 0116583 
23 ei41j{26 
26 0159534 
28 0171806 
30 0184078 
31 0190214 
u2 2257709 
58 0355884 

I BIGMA(I) 
68 044720277 
73 646019484 
76 046633016 
17 e47246608 

NSALFRa 

DEPTH OF WATER AT GAUGE SITEm 

AVGy« 222246 AvVG2a 

SUM)s 0299 5UM2g 

WSUMI a9 2084 WSUMBa 

RATIO{# 30579 RATIU2e 

SUMENE 031019470 


BREAKING WAyE HEIGHT Hb 
BREAKING WAv& CELERITY CBs 


RSODDe 290048 RSHFRQa 
PLPQ8e 125.7135 PLNEGa 
PLNETe 109.0401 

Figure 8. 


Islands Harbor.-—Continued 
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TIME 
1800 
1800 


FMD8O(1) 
0000169 
2000099 
2000007 
2000013 
©000066 
0000164 
2000038 
2000000 
2000168 
2000201 
0000137 
0000114 
2000055 
0000061 
2000004 
2000072 
2000008 
0000020 
2000006 
20000288 
0000064 


PCT 
002575229 
0048453916 
cOS14i7S3 
004237300 


203 
24.0 
20.808 

0293 
2078 
3.741 


3.61 
12.22 


THETACT) 
075283813 
065863R27 
069562037 
o7309U726 


3742 RSLFRUa 00114 


25.6734 


Three examples of output for wave gage pair at Channel 


Example 3 


GAUGE NO, MONTH DAY TIME 
Sit u 26 2000 
312 7 26 2000 

if SIGMA(1) PMpSaQ(y) 
1 2006136 2000930 
2 vOl2272 2000397 
4 0018408 2000179 
u 2024544 2000082 
b 0042951 0000013 
8 0049087 2000081 
10 0061359 2000084 
1a 0073031 0000074 
13 eO7T9767 2000085 
17 s/O43t1 - 2000080 
19 0116563 0000008 
23 o141126 00000K1 
32 0196350 2000009 
34 0806621 0000100 
3§ 08144797 2000104 
36 02208693 0000019 
i } 0aS7709 2000044 
53 o337476 0000136 
11 0435651 0001082 
1 SIGMaA(]) PCT THEeTa(t) 

NGALFRe ao02 

DEPTH OF WATER AT GAUGE SITEs 23.8 

AVG,a a, .702 AVG2a 20.487 

8UM;s 0253 8UMa@ 0266 

WEUMi a 0073 WEUMae 2081 

RATIVIs 30460 RATIV2a 3,299 

SUMEN® 033101957 

BREAKING WAVE HEIGHT Hes 3,66 

BREAKING WAyE CELBRITy CBe 12.29 

R80DDa 0095 RBHFRUe 0867S ROLFRUe 0142 

PLPO8s 135.5367 PLNEGR 040.2297 

PLNETs 95.3070 


Figure 8. Three examples of output for wave gage pair at Channel 
Islands Harbor.——Continued 
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